/*
 * CUDA_JacobiRelax.h
 *
 *  Created on: Mar 22, 2011
 *      Author: curt
 */

#ifndef CUDA_JACOBIRELAX_H_
#define CUDA_JACOBIRELAX_H_

#include "CUDA_vec.h"

class CUDA_JacobiRelax {
public:
	typedef CUDA_vec<float> Vector;

	CUDA_JacobiRelax(int num_iterations, float relax_param):num_iter(num_iterations),omega(relax_param){}

	template<typename MatrixType>
	void operator()(const MatrixType & A, Vector & x, const Vector & b)
	{
		CUDA_vec<float> invd(b.getSize());
		CUDA_vec<float> y(b.getSize());
		A.getDiagonal(invd);
		invd.reciprocal();

		for(int i=0; i<num_iter; i++)
		{
			//y = (b - A*x) ./ diagonal
			A(x, y);
			y.AXPBY(-1, 1, b);
			y*= invd;

			//x = x + omega * y
			x.AXPBY(1,omega, y);
		}

	}

private:
	int num_iter;
	float omega;
};


#endif /* CUDA_JACOBIRELAX_H_ */
